home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dtrsyl.f < prev    next >
Text File  |  1996-07-19  |  34KB  |  915 lines
  1.       SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
  2.      $                   LDC, SCALE, INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     March 31, 1993
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          TRANA, TRANB
  11.       INTEGER            INFO, ISGN, LDA, LDB, LDC, M, N
  12.       DOUBLE PRECISION   SCALE
  13. *     ..
  14. *     .. Array Arguments ..
  15.       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  DTRSYL solves the real Sylvester matrix equation:
  22. *
  23. *     op(A)*X + X*op(B) = scale*C or
  24. *     op(A)*X - X*op(B) = scale*C,
  25. *
  26. *  where op(A) = A or A**T, and  A and B are both upper quasi-
  27. *  triangular. A is M-by-M and B is N-by-N; the right hand side C and
  28. *  the solution X are M-by-N; and scale is an output scale factor, set
  29. *  <= 1 to avoid overflow in X.
  30. *
  31. *  A and B must be in Schur canonical form (as returned by DHSEQR), that
  32. *  is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
  33. *  each 2-by-2 diagonal block has its diagonal elements equal and its
  34. *  off-diagonal elements of opposite sign.
  35. *
  36. *  Arguments
  37. *  =========
  38. *
  39. *  TRANA   (input) CHARACTER*1
  40. *          Specifies the option op(A):
  41. *          = 'N': op(A) = A    (No transpose)
  42. *          = 'T': op(A) = A**T (Transpose)
  43. *          = 'C': op(A) = A**H (Conjugate transpose = Transpose)
  44. *
  45. *  TRANB   (input) CHARACTER*1
  46. *          Specifies the option op(B):
  47. *          = 'N': op(B) = B    (No transpose)
  48. *          = 'T': op(B) = B**T (Transpose)
  49. *          = 'C': op(B) = B**H (Conjugate transpose = Transpose)
  50. *
  51. *  ISGN    (input) INTEGER
  52. *          Specifies the sign in the equation:
  53. *          = +1: solve op(A)*X + X*op(B) = scale*C
  54. *          = -1: solve op(A)*X - X*op(B) = scale*C
  55. *
  56. *  M       (input) INTEGER
  57. *          The order of the matrix A, and the number of rows in the
  58. *          matrices X and C. M >= 0.
  59. *
  60. *  N       (input) INTEGER
  61. *          The order of the matrix B, and the number of columns in the
  62. *          matrices X and C. N >= 0.
  63. *
  64. *  A       (input) DOUBLE PRECISION array, dimension (LDA,M)
  65. *          The upper quasi-triangular matrix A, in Schur canonical form.
  66. *
  67. *  LDA     (input) INTEGER
  68. *          The leading dimension of the array A. LDA >= max(1,M).
  69. *
  70. *  B       (input) DOUBLE PRECISION array, dimension (LDB,N)
  71. *          The upper quasi-triangular matrix B, in Schur canonical form.
  72. *
  73. *  LDB     (input) INTEGER
  74. *          The leading dimension of the array B. LDB >= max(1,N).
  75. *
  76. *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
  77. *          On entry, the M-by-N right hand side matrix C.
  78. *          On exit, C is overwritten by the solution matrix X.
  79. *
  80. *  LDC     (input) INTEGER
  81. *          The leading dimension of the array C. LDC >= max(1,M)
  82. *
  83. *  SCALE   (output) DOUBLE PRECISION
  84. *          The scale factor, scale, set <= 1 to avoid overflow in X.
  85. *
  86. *  INFO    (output) INTEGER
  87. *          = 0: successful exit
  88. *          < 0: if INFO = -i, the i-th argument had an illegal value
  89. *          = 1: A and B have common or very close eigenvalues; perturbed
  90. *               values were used to solve the equation (but the matrices
  91. *               A and B are unchanged).
  92. *
  93. *  =====================================================================
  94. *
  95. *     .. Parameters ..
  96.       DOUBLE PRECISION   ZERO, ONE
  97.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  98. *     ..
  99. *     .. Local Scalars ..
  100.       LOGICAL            NOTRNA, NOTRNB
  101.       INTEGER            IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT
  102.       DOUBLE PRECISION   A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
  103.      $                   SMLNUM, SUML, SUMR, XNORM
  104. *     ..
  105. *     .. Local Arrays ..
  106.       DOUBLE PRECISION   DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
  107. *     ..
  108. *     .. External Functions ..
  109.       LOGICAL            LSAME
  110.       DOUBLE PRECISION   DDOT, DLAMCH, DLANGE
  111.       EXTERNAL           LSAME, DDOT, DLAMCH, DLANGE
  112. *     ..
  113. *     .. External Subroutines ..
  114.       EXTERNAL           DLABAD, DLALN2, DLASY2, DSCAL, XERBLA
  115. *     ..
  116. *     .. Intrinsic Functions ..
  117.       INTRINSIC          ABS, DBLE, MAX, MIN
  118. *     ..
  119. *     .. Executable Statements ..
  120. *
  121. *     Decode and Test input parameters
  122. *
  123.       NOTRNA = LSAME( TRANA, 'N' )
  124.       NOTRNB = LSAME( TRANB, 'N' )
  125. *
  126.       INFO = 0
  127.       IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT.
  128.      $    LSAME( TRANA, 'C' ) ) THEN
  129.          INFO = -1
  130.       ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT.
  131.      $         LSAME( TRANB, 'C' ) ) THEN
  132.          INFO = -2
  133.       ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
  134.          INFO = -3
  135.       ELSE IF( M.LT.0 ) THEN
  136.          INFO = -4
  137.       ELSE IF( N.LT.0 ) THEN
  138.          INFO = -5
  139.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  140.          INFO = -7
  141.       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  142.          INFO = -9
  143.       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
  144.          INFO = -11
  145.       END IF
  146.       IF( INFO.NE.0 ) THEN
  147.          CALL XERBLA( 'DTRSYL', -INFO )
  148.          RETURN
  149.       END IF
  150. *
  151. *     Quick return if possible
  152. *
  153.       IF( M.EQ.0 .OR. N.EQ.0 )
  154.      $   RETURN
  155. *
  156. *     Set constants to control overflow
  157. *
  158.       EPS = DLAMCH( 'P' )
  159.       SMLNUM = DLAMCH( 'S' )
  160.       BIGNUM = ONE / SMLNUM
  161.       CALL DLABAD( SMLNUM, BIGNUM )
  162.       SMLNUM = SMLNUM*DBLE( M*N ) / EPS
  163.       BIGNUM = ONE / SMLNUM
  164. *
  165.       SMIN = MAX( SMLNUM, EPS*DLANGE( 'M', M, M, A, LDA, DUM ),
  166.      $       EPS*DLANGE( 'M', N, N, B, LDB, DUM ) )
  167. *
  168.       SCALE = ONE
  169.       SGN = ISGN
  170. *
  171.       IF( NOTRNA .AND. NOTRNB ) THEN
  172. *
  173. *        Solve    A*X + ISGN*X*B = scale*C.
  174. *
  175. *        The (K,L)th block of X is determined starting from
  176. *        bottom-left corner column by column by
  177. *
  178. *         A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
  179. *
  180. *        Where
  181. *                  M                         L-1
  182. *        R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
  183. *                I=K+1                       J=1
  184. *
  185. *        Start column loop (index = L)
  186. *        L1 (L2) : column index of the first (first) row of X(K,L).
  187. *
  188.          LNEXT = 1
  189.          DO 60 L = 1, N
  190.             IF( L.LT.LNEXT )
  191.      $         GO TO 60
  192.             IF( L.EQ.N ) THEN
  193.                L1 = L
  194.                L2 = L
  195.             ELSE
  196.                IF( B( L+1, L ).NE.ZERO ) THEN
  197.                   L1 = L
  198.                   L2 = L + 1
  199.                   LNEXT = L + 2
  200.                ELSE
  201.                   L1 = L
  202.                   L2 = L
  203.                   LNEXT = L + 1
  204.                END IF
  205.             END IF
  206. *
  207. *           Start row loop (index = K)
  208. *           K1 (K2): row index of the first (last) row of X(K,L).
  209. *
  210.             KNEXT = M
  211.             DO 50 K = M, 1, -1
  212.                IF( K.GT.KNEXT )
  213.      $            GO TO 50
  214.                IF( K.EQ.1 ) THEN
  215.                   K1 = K
  216.                   K2 = K
  217.                ELSE
  218.                   IF( A( K, K-1 ).NE.ZERO ) THEN
  219.                      K1 = K - 1
  220.                      K2 = K
  221.                      KNEXT = K - 2
  222.                   ELSE
  223.                      K1 = K
  224.                      K2 = K
  225.                      KNEXT = K - 1
  226.                   END IF
  227.                END IF
  228. *
  229.                IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
  230.                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  231.      $                   C( MIN( K1+1, M ), L1 ), 1 )
  232.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  233.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  234.                   SCALOC = ONE
  235. *
  236.                   A11 = A( K1, K1 ) + SGN*B( L1, L1 )
  237.                   DA11 = ABS( A11 )
  238.                   IF( DA11.LE.SMIN ) THEN
  239.                      A11 = SMIN
  240.                      DA11 = SMIN
  241.                      INFO = 1
  242.                   END IF
  243.                   DB = ABS( VEC( 1, 1 ) )
  244.                   IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  245.                      IF( DB.GT.BIGNUM*DA11 )
  246.      $                  SCALOC = ONE / DB
  247.                   END IF
  248.                   X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
  249. *
  250.                   IF( SCALOC.NE.ONE ) THEN
  251.                      DO 10 J = 1, N
  252.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  253.    10                CONTINUE
  254.                      SCALE = SCALE*SCALOC
  255.                   END IF
  256.                   C( K1, L1 ) = X( 1, 1 )
  257. *
  258.                ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
  259. *
  260.                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  261.      $                   C( MIN( K2+1, M ), L1 ), 1 )
  262.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  263.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  264. *
  265.                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  266.      $                   C( MIN( K2+1, M ), L1 ), 1 )
  267.                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
  268.                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  269. *
  270.                   CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
  271.      $                         LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
  272.      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
  273.                   IF( IERR.NE.0 )
  274.      $               INFO = 1
  275. *
  276.                   IF( SCALOC.NE.ONE ) THEN
  277.                      DO 20 J = 1, N
  278.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  279.    20                CONTINUE
  280.                      SCALE = SCALE*SCALOC
  281.                   END IF
  282.                   C( K1, L1 ) = X( 1, 1 )
  283.                   C( K2, L1 ) = X( 2, 1 )
  284. *
  285.                ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
  286. *
  287.                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  288.      $                   C( MIN( K1+1, M ), L1 ), 1 )
  289.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  290.                   VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
  291. *
  292.                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  293.      $                   C( MIN( K1+1, M ), L2 ), 1 )
  294.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
  295.                   VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
  296. *
  297.                   CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ),
  298.      $                         LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
  299.      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
  300.                   IF( IERR.NE.0 )
  301.      $               INFO = 1
  302. *
  303.                   IF( SCALOC.NE.ONE ) THEN
  304.                      DO 30 J = 1, N
  305.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  306.    30                CONTINUE
  307.                      SCALE = SCALE*SCALOC
  308.                   END IF
  309.                   C( K1, L1 ) = X( 1, 1 )
  310.                   C( K1, L2 ) = X( 2, 1 )
  311. *
  312.                ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
  313. *
  314.                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  315.      $                   C( MIN( K2+1, M ), L1 ), 1 )
  316.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  317.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  318. *
  319.                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  320.      $                   C( MIN( K2+1, M ), L2 ), 1 )
  321.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
  322.                   VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
  323. *
  324.                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  325.      $                   C( MIN( K2+1, M ), L1 ), 1 )
  326.                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
  327.                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  328. *
  329.                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  330.      $                   C( MIN( K2+1, M ), L2 ), 1 )
  331.                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 )
  332.                   VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
  333. *
  334.                   CALL DLASY2( .FALSE., .FALSE., ISGN, 2, 2,
  335.      $                         A( K1, K1 ), LDA, B( L1, L1 ), LDB, VEC,
  336.      $                         2, SCALOC, X, 2, XNORM, IERR )
  337.                   IF( IERR.NE.0 )
  338.      $               INFO = 1
  339. *
  340.                   IF( SCALOC.NE.ONE ) THEN
  341.                      DO 40 J = 1, N
  342.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  343.    40                CONTINUE
  344.                      SCALE = SCALE*SCALOC
  345.                   END IF
  346.                   C( K1, L1 ) = X( 1, 1 )
  347.                   C( K1, L2 ) = X( 1, 2 )
  348.                   C( K2, L1 ) = X( 2, 1 )
  349.                   C( K2, L2 ) = X( 2, 2 )
  350.                END IF
  351. *
  352.    50       CONTINUE
  353. *
  354.    60    CONTINUE
  355. *
  356.       ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
  357. *
  358. *        Solve    A' *X + ISGN*X*B = scale*C.
  359. *
  360. *        The (K,L)th block of X is determined starting from
  361. *        upper-left corner column by column by
  362. *
  363. *          A(K,K)'*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
  364. *
  365. *        Where
  366. *                   K-1                        L-1
  367. *          R(K,L) = SUM [A(I,K)'*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
  368. *                   I=1                        J=1
  369. *
  370. *        Start column loop (index = L)
  371. *        L1 (L2): column index of the first (last) row of X(K,L)
  372. *
  373.          LNEXT = 1
  374.          DO 120 L = 1, N
  375.             IF( L.LT.LNEXT )
  376.      $         GO TO 120
  377.             IF( L.EQ.N ) THEN
  378.                L1 = L
  379.                L2 = L
  380.             ELSE
  381.                IF( B( L+1, L ).NE.ZERO ) THEN
  382.                   L1 = L
  383.                   L2 = L + 1
  384.                   LNEXT = L + 2
  385.                ELSE
  386.                   L1 = L
  387.                   L2 = L
  388.                   LNEXT = L + 1
  389.                END IF
  390.             END IF
  391. *
  392. *           Start row loop (index = K)
  393. *           K1 (K2): row index of the first (last) row of X(K,L)
  394. *
  395.             KNEXT = 1
  396.             DO 110 K = 1, M
  397.                IF( K.LT.KNEXT )
  398.      $            GO TO 110
  399.                IF( K.EQ.M ) THEN
  400.                   K1 = K
  401.                   K2 = K
  402.                ELSE
  403.                   IF( A( K+1, K ).NE.ZERO ) THEN
  404.                      K1 = K
  405.                      K2 = K + 1
  406.                      KNEXT = K + 2
  407.                   ELSE
  408.                      K1 = K
  409.                      K2 = K
  410.                      KNEXT = K + 1
  411.                   END IF
  412.                END IF
  413. *
  414.                IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
  415.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  416.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  417.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  418.                   SCALOC = ONE
  419. *
  420.                   A11 = A( K1, K1 ) + SGN*B( L1, L1 )
  421.                   DA11 = ABS( A11 )
  422.                   IF( DA11.LE.SMIN ) THEN
  423.                      A11 = SMIN
  424.                      DA11 = SMIN
  425.                      INFO = 1
  426.                   END IF
  427.                   DB = ABS( VEC( 1, 1 ) )
  428.                   IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  429.                      IF( DB.GT.BIGNUM*DA11 )
  430.      $                  SCALOC = ONE / DB
  431.                   END IF
  432.                   X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
  433. *
  434.                   IF( SCALOC.NE.ONE ) THEN
  435.                      DO 70 J = 1, N
  436.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  437.    70                CONTINUE
  438.                      SCALE = SCALE*SCALOC
  439.                   END IF
  440.                   C( K1, L1 ) = X( 1, 1 )
  441. *
  442.                ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
  443. *
  444.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  445.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  446.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  447. *
  448.                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
  449.                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
  450.                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  451. *
  452.                   CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
  453.      $                         LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
  454.      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
  455.                   IF( IERR.NE.0 )
  456.      $               INFO = 1
  457. *
  458.                   IF( SCALOC.NE.ONE ) THEN
  459.                      DO 80 J = 1, N
  460.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  461.    80                CONTINUE
  462.                      SCALE = SCALE*SCALOC
  463.                   END IF
  464.                   C( K1, L1 ) = X( 1, 1 )
  465.                   C( K2, L1 ) = X( 2, 1 )
  466. *
  467.                ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
  468. *
  469.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  470.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  471.                   VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
  472. *
  473.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
  474.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
  475.                   VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
  476. *
  477.                   CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ),
  478.      $                         LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
  479.      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
  480.                   IF( IERR.NE.0 )
  481.      $               INFO = 1
  482. *
  483.                   IF( SCALOC.NE.ONE ) THEN
  484.                      DO 90 J = 1, N
  485.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  486.    90                CONTINUE
  487.                      SCALE = SCALE*SCALOC
  488.                   END IF
  489.                   C( K1, L1 ) = X( 1, 1 )
  490.                   C( K1, L2 ) = X( 2, 1 )
  491. *
  492.                ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
  493. *
  494.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  495.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  496.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  497. *
  498.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
  499.                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
  500.                   VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
  501. *
  502.                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
  503.                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
  504.                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  505. *
  506.                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 )
  507.                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 )
  508.                   VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
  509. *
  510.                   CALL DLASY2( .TRUE., .FALSE., ISGN, 2, 2, A( K1, K1 ),
  511.      $                         LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
  512.      $                         2, XNORM, IERR )
  513.                   IF( IERR.NE.0 )
  514.      $               INFO = 1
  515. *
  516.                   IF( SCALOC.NE.ONE ) THEN
  517.                      DO 100 J = 1, N
  518.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  519.   100                CONTINUE
  520.                      SCALE = SCALE*SCALOC
  521.                   END IF
  522.                   C( K1, L1 ) = X( 1, 1 )
  523.                   C( K1, L2 ) = X( 1, 2 )
  524.                   C( K2, L1 ) = X( 2, 1 )
  525.                   C( K2, L2 ) = X( 2, 2 )
  526.                END IF
  527. *
  528.   110       CONTINUE
  529.   120    CONTINUE
  530. *
  531.       ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
  532. *
  533. *        Solve    A'*X + ISGN*X*B' = scale*C.
  534. *
  535. *        The (K,L)th block of X is determined starting from
  536. *        top-right corner column by column by
  537. *
  538. *           A(K,K)'*X(K,L) + ISGN*X(K,L)*B(L,L)' = C(K,L) - R(K,L)
  539. *
  540. *        Where
  541. *                     K-1                          N
  542. *            R(K,L) = SUM [A(I,K)'*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)'].
  543. *                     I=1                        J=L+1
  544. *
  545. *        Start column loop (index = L)
  546. *        L1 (L2): column index of the first (last) row of X(K,L)
  547. *
  548.          LNEXT = N
  549.          DO 180 L = N, 1, -1
  550.             IF( L.GT.LNEXT )
  551.      $         GO TO 180
  552.             IF( L.EQ.1 ) THEN
  553.                L1 = L
  554.                L2 = L
  555.             ELSE
  556.                IF( B( L, L-1 ).NE.ZERO ) THEN
  557.                   L1 = L - 1
  558.                   L2 = L
  559.                   LNEXT = L - 2
  560.                ELSE
  561.                   L1 = L
  562.                   L2 = L
  563.                   LNEXT = L - 1
  564.                END IF
  565.             END IF
  566. *
  567. *           Start row loop (index = K)
  568. *           K1 (K2): row index of the first (last) row of X(K,L)
  569. *
  570.             KNEXT = 1
  571.             DO 170 K = 1, M
  572.                IF( K.LT.KNEXT )
  573.      $            GO TO 170
  574.                IF( K.EQ.M ) THEN
  575.                   K1 = K
  576.                   K2 = K
  577.                ELSE
  578.                   IF( A( K+1, K ).NE.ZERO ) THEN
  579.                      K1 = K
  580.                      K2 = K + 1
  581.                      KNEXT = K + 2
  582.                   ELSE
  583.                      K1 = K
  584.                      K2 = K
  585.                      KNEXT = K + 1
  586.                   END IF
  587.                END IF
  588. *
  589.                IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
  590.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  591.                   SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC,
  592.      $                   B( L1, MIN( L1+1, N ) ), LDB )
  593.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  594.                   SCALOC = ONE
  595. *
  596.                   A11 = A( K1, K1 ) + SGN*B( L1, L1 )
  597.                   DA11 = ABS( A11 )
  598.                   IF( DA11.LE.SMIN ) THEN
  599.                      A11 = SMIN
  600.                      DA11 = SMIN
  601.                      INFO = 1
  602.                   END IF
  603.                   DB = ABS( VEC( 1, 1 ) )
  604.                   IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  605.                      IF( DB.GT.BIGNUM*DA11 )
  606.      $                  SCALOC = ONE / DB
  607.                   END IF
  608.                   X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
  609. *
  610.                   IF( SCALOC.NE.ONE ) THEN
  611.                      DO 130 J = 1, N
  612.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  613.   130                CONTINUE
  614.                      SCALE = SCALE*SCALOC
  615.                   END IF
  616.                   C( K1, L1 ) = X( 1, 1 )
  617. *
  618.                ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
  619. *
  620.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  621.                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  622.      $                   B( L1, MIN( L2+1, N ) ), LDB )
  623.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  624. *
  625.                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
  626.                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  627.      $                   B( L1, MIN( L2+1, N ) ), LDB )
  628.                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  629. *
  630.                   CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
  631.      $                         LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
  632.      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
  633.                   IF( IERR.NE.0 )
  634.      $               INFO = 1
  635. *
  636.                   IF( SCALOC.NE.ONE ) THEN
  637.                      DO 140 J = 1, N
  638.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  639.   140                CONTINUE
  640.                      SCALE = SCALE*SCALOC
  641.                   END IF
  642.                   C( K1, L1 ) = X( 1, 1 )
  643.                   C( K2, L1 ) = X( 2, 1 )
  644. *
  645.                ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
  646. *
  647.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  648.                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  649.      $                   B( L1, MIN( L2+1, N ) ), LDB )
  650.                   VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
  651. *
  652.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
  653.                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  654.      $                   B( L2, MIN( L2+1, N ) ), LDB )
  655.                   VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
  656. *
  657.                   CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ),
  658.      $                         LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
  659.      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
  660.                   IF( IERR.NE.0 )
  661.      $               INFO = 1
  662. *
  663.                   IF( SCALOC.NE.ONE ) THEN
  664.                      DO 150 J = 1, N
  665.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  666.   150                CONTINUE
  667.                      SCALE = SCALE*SCALOC
  668.                   END IF
  669.                   C( K1, L1 ) = X( 1, 1 )
  670.                   C( K1, L2 ) = X( 2, 1 )
  671. *
  672.                ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
  673. *
  674.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  675.                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  676.      $                   B( L1, MIN( L2+1, N ) ), LDB )
  677.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  678. *
  679.                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
  680.                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  681.      $                   B( L2, MIN( L2+1, N ) ), LDB )
  682.                   VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
  683. *
  684.                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
  685.                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  686.      $                   B( L1, MIN( L2+1, N ) ), LDB )
  687.                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  688. *
  689.                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 )
  690.                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  691.      $                   B( L2, MIN( L2+1, N ) ), LDB )
  692.                   VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
  693. *
  694.                   CALL DLASY2( .TRUE., .TRUE., ISGN, 2, 2, A( K1, K1 ),
  695.      $                         LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
  696.      $                         2, XNORM, IERR )
  697.                   IF( IERR.NE.0 )
  698.      $               INFO = 1
  699. *
  700.                   IF( SCALOC.NE.ONE ) THEN
  701.                      DO 160 J = 1, N
  702.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  703.   160                CONTINUE
  704.                      SCALE = SCALE*SCALOC
  705.                   END IF
  706.                   C( K1, L1 ) = X( 1, 1 )
  707.                   C( K1, L2 ) = X( 1, 2 )
  708.                   C( K2, L1 ) = X( 2, 1 )
  709.                   C( K2, L2 ) = X( 2, 2 )
  710.                END IF
  711. *
  712.   170       CONTINUE
  713.   180    CONTINUE
  714. *
  715.       ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
  716. *
  717. *        Solve    A*X + ISGN*X*B' = scale*C.
  718. *
  719. *        The (K,L)th block of X is determined starting from
  720. *        bottom-right corner column by column by
  721. *
  722. *            A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)' = C(K,L) - R(K,L)
  723. *
  724. *        Where
  725. *                      M                          N
  726. *            R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)'].
  727. *                    I=K+1                      J=L+1
  728. *
  729. *        Start column loop (index = L)
  730. *        L1 (L2): column index of the first (last) row of X(K,L)
  731. *
  732.          LNEXT = N
  733.          DO 240 L = N, 1, -1
  734.             IF( L.GT.LNEXT )
  735.      $         GO TO 240
  736.             IF( L.EQ.1 ) THEN
  737.                L1 = L
  738.                L2 = L
  739.             ELSE
  740.                IF( B( L, L-1 ).NE.ZERO ) THEN
  741.                   L1 = L - 1
  742.                   L2 = L
  743.                   LNEXT = L - 2
  744.                ELSE
  745.                   L1 = L
  746.                   L2 = L
  747.                   LNEXT = L - 1
  748.                END IF
  749.             END IF
  750. *
  751. *           Start row loop (index = K)
  752. *           K1 (K2): row index of the first (last) row of X(K,L)
  753. *
  754.             KNEXT = M
  755.             DO 230 K = M, 1, -1
  756.                IF( K.GT.KNEXT )
  757.      $            GO TO 230
  758.                IF( K.EQ.1 ) THEN
  759.                   K1 = K
  760.                   K2 = K
  761.                ELSE
  762.                   IF( A( K, K-1 ).NE.ZERO ) THEN
  763.                      K1 = K - 1
  764.                      K2 = K
  765.                      KNEXT = K - 2
  766.                   ELSE
  767.                      K1 = K
  768.                      K2 = K
  769.                      KNEXT = K - 1
  770.                   END IF
  771.                END IF
  772. *
  773.                IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
  774.                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  775.      $                   C( MIN( K1+1, M ), L1 ), 1 )
  776.                   SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC,
  777.      $                   B( L1, MIN( L1+1, N ) ), LDB )
  778.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  779.                   SCALOC = ONE
  780. *
  781.                   A11 = A( K1, K1 ) + SGN*B( L1, L1 )
  782.                   DA11 = ABS( A11 )
  783.                   IF( DA11.LE.SMIN ) THEN
  784.                      A11 = SMIN
  785.                      DA11 = SMIN
  786.                      INFO = 1
  787.                   END IF
  788.                   DB = ABS( VEC( 1, 1 ) )
  789.                   IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  790.                      IF( DB.GT.BIGNUM*DA11 )
  791.      $                  SCALOC = ONE / DB
  792.                   END IF
  793.                   X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
  794. *
  795.                   IF( SCALOC.NE.ONE ) THEN
  796.                      DO 190 J = 1, N
  797.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  798.   190                CONTINUE
  799.                      SCALE = SCALE*SCALOC
  800.                   END IF
  801.                   C( K1, L1 ) = X( 1, 1 )
  802. *
  803.                ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
  804. *
  805.                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  806.      $                   C( MIN( K2+1, M ), L1 ), 1 )
  807.                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  808.      $                   B( L1, MIN( L2+1, N ) ), LDB )
  809.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  810. *
  811.                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  812.      $                   C( MIN( K2+1, M ), L1 ), 1 )
  813.                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  814.      $                   B( L1, MIN( L2+1, N ) ), LDB )
  815.                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  816. *
  817.                   CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
  818.      $                         LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
  819.      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
  820.                   IF( IERR.NE.0 )
  821.      $               INFO = 1
  822. *
  823.                   IF( SCALOC.NE.ONE ) THEN
  824.                      DO 200 J = 1, N
  825.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  826.   200                CONTINUE
  827.                      SCALE = SCALE*SCALOC
  828.                   END IF
  829.                   C( K1, L1 ) = X( 1, 1 )
  830.                   C( K2, L1 ) = X( 2, 1 )
  831. *
  832.                ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
  833. *
  834.                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  835.      $                   C( MIN( K1+1, M ), L1 ), 1 )
  836.                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  837.      $                   B( L1, MIN( L2+1, N ) ), LDB )
  838.                   VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
  839. *
  840.                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  841.      $                   C( MIN( K1+1, M ), L2 ), 1 )
  842.                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  843.      $                   B( L2, MIN( L2+1, N ) ), LDB )
  844.                   VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
  845. *
  846.                   CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ),
  847.      $                         LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
  848.      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
  849.                   IF( IERR.NE.0 )
  850.      $               INFO = 1
  851. *
  852.                   IF( SCALOC.NE.ONE ) THEN
  853.                      DO 210 J = 1, N
  854.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  855.   210                CONTINUE
  856.                      SCALE = SCALE*SCALOC
  857.                   END IF
  858.                   C( K1, L1 ) = X( 1, 1 )
  859.                   C( K1, L2 ) = X( 2, 1 )
  860. *
  861.                ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
  862. *
  863.                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  864.      $                   C( MIN( K2+1, M ), L1 ), 1 )
  865.                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  866.      $                   B( L1, MIN( L2+1, N ) ), LDB )
  867.                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  868. *
  869.                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  870.      $                   C( MIN( K2+1, M ), L2 ), 1 )
  871.                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  872.      $                   B( L2, MIN( L2+1, N ) ), LDB )
  873.                   VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
  874. *
  875.                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  876.      $                   C( MIN( K2+1, M ), L1 ), 1 )
  877.                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  878.      $                   B( L1, MIN( L2+1, N ) ), LDB )
  879.                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  880. *
  881.                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  882.      $                   C( MIN( K2+1, M ), L2 ), 1 )
  883.                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  884.      $                   B( L2, MIN( L2+1, N ) ), LDB )
  885.                   VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
  886. *
  887.                   CALL DLASY2( .FALSE., .TRUE., ISGN, 2, 2, A( K1, K1 ),
  888.      $                         LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
  889.      $                         2, XNORM, IERR )
  890.                   IF( IERR.NE.0 )
  891.      $               INFO = 1
  892. *
  893.                   IF( SCALOC.NE.ONE ) THEN
  894.                      DO 220 J = 1, N
  895.                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  896.   220                CONTINUE
  897.                      SCALE = SCALE*SCALOC
  898.                   END IF
  899.                   C( K1, L1 ) = X( 1, 1 )
  900.                   C( K1, L2 ) = X( 1, 2 )
  901.                   C( K2, L1 ) = X( 2, 1 )
  902.                   C( K2, L2 ) = X( 2, 2 )
  903.                END IF
  904. *
  905.   230       CONTINUE
  906.   240    CONTINUE
  907. *
  908.       END IF
  909. *
  910.       RETURN
  911. *
  912. *     End of DTRSYL
  913. *
  914.       END
  915.